Simulators are an important tool for the design and control of robots, vehicles, and other cyber-physical systems. Many simulators are based on multibody system dynamics in maximal coordinates with joints and contacts represented as kinematic constraints and complementarity conditions. This allows for fast simulation using time-implicit integration but involves model and simulation parameters that are not derivable from first-principle models or identifiable by direct measurements. If these parameters are not accurately determined, the result of simulations transfer poorly to the real world. Therefore, there is a need for scalable methods and highly automated tools for the task of parameter identification. We present a method for offline parameter estimation of multibody dynamics simulators utilizing maximal coordinates, constraint stabilization and regularization [2]. This formulation is essential to handle frictional contacts, kinematic loops and non-ideal constraints, but poses a problem for observability – joint angles no longer map uniquely to body poses. To this end, we solve the simultaneous state and parameter estimation problem based on a non-linear least squares formulation, motivated by maximum likelihood estimation, and similar to filtering methods. Derivatives are computed using forward-mode automatic differentiation, with special treatment of the linear complementarity problem and rotations. The feasibility is demonstrated through the identification of parameters of a Furuta pendulum, an underactuated mechatronic system composed of three rigid bodies, three joints, and an electric motor. The parameters subject for identification include center of mass, inertia, dry friction, viscous friction, direction of gravity, motor gain, joint clearance, and joint compliance. Results show high agreement when the measurements include scenarios where the parameters are properly excited. We discuss the challenges and possibilities with automating the process and extending the method to include non-smooth dynamics.